library(rethinking)
Loading required package: rstan
Loading required package: StanHeaders
Loading required package: ggplot2
rstan (Version 2.21.2, GitRev: 2e1f913d3ca3)
For execution on a local, multicore CPU with excess RAM we recommend calling
options(mc.cores = parallel::detectCores()).
To avoid recompilation of unchanged Stan programs, we recommend calling
rstan_options(auto_write = TRUE)
Loading required package: parallel
rethinking (Version 2.11)
Attaching package: ‘rethinking’
The following object is masked from ‘package:stats’:
rstudent
library(brms)
Loading required package: Rcpp
Registered S3 method overwritten by 'htmlwidgets':
method from
print.htmlwidget tools:rstudio
Loading 'brms' package (version 2.13.5). Useful instructions
can be found by typing help('brms'). A more detailed introduction
to the package is available through vignette('brms_overview').
Attaching package: ‘brms’
The following objects are masked from ‘package:rethinking’:
LOO, stancode, WAIC
The following object is masked from ‘package:rstan’:
loo
The following object is masked from ‘package:stats’:
ar
library(tidyverse)
Registered S3 methods overwritten by 'dbplyr':
method from
print.tbl_lazy
print.tbl_sql
── Attaching packages ───────────────────────────────────────────────────────────────── tidyverse 1.3.0 ──
✓ tibble 3.0.3 ✓ dplyr 1.0.0
✓ tidyr 1.1.0 ✓ stringr 1.4.0
✓ readr 1.3.1 ✓ forcats 0.5.0
✓ purrr 0.3.4
── Conflicts ──────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
x tidyr::extract() masks rstan::extract()
x dplyr::filter() masks stats::filter()
x dplyr::lag() masks stats::lag()
x purrr::map() masks rethinking::map()
germ <- read_csv("../Dimensions/hydrothermal-all-species/data/light_round1_tall.csv") %>% filter(wps==0) %>%
Warning message:
Unknown or uninitialised column: `germ`.
select(pops, temps, total_seeds, germ, day, cumulative_germ)
Parsed with column specification:
cols(
pops = col_character(),
temps = col_double(),
wps = col_double(),
date = col_character(),
total_seeds = col_double(),
germ = col_double(),
start_date = col_date(format = ""),
census_date = col_date(format = ""),
day = col_double(),
cumulative_germ = col_double(),
cumulative_prop_germ = col_double()
)
germ
what if need one event per row:
one_per_row <- function(df) {
total_seed <- max(df$total_seeds, sum(df$germ))
newdata <- tibble(id=1:total_seed, germ=0, day=max(df$day))
df <- df %>% filter(germ>0)
count <- 1
if (nrow(df) > 0) {
for (i in 1:nrow(df)) { # we look at each row of the df where germination occured
for (j in 1:df$germ[i]) { # now update the newdata to reflect the germiantion of each seed
newdata$germ[count] <- 1
newdata$day[count]=df$day[i]
count <- count+1 # count keeps track of which individual we are at in the new data
} # for j
} # for i
} # if
return(newdata)
}
germone <- germ %>% group_by(pops, temps) %>%
select(-cumulative_germ) %>% # not needed in this encoding (I think...in any case would need to be recalculated)
nest() %>%
mutate(newdata=map(data, one_per_row)) %>%
select(-data) %>%
unnest(newdata)
germone
##1
Is effect of temp ~linear on cum germ?
germ %>% filter(pops=="STDI", day==28) %>%
Warning message:
Unknown or uninitialised column: `germ`.
ggplot(aes(x=temps,y=cumulative_germ)) +
geom_col()
no, so we can’t treat temp as a continuous linear predictor
germ.stdi <- germone %>% filter(pops=="STDI") %>% select(-pops)
Adding missing grouping variables: `pops`
germ.stdi
Try a censored lambda model
d <- list(germ=germ.stdi$germ, temps=as.numeric(as.factor(germ.stdi$temps)), day=germ.stdi$day)
Warning messages:
1: Unknown or uninitialised column: `germ`.
2: Unknown or uninitialised column: `germ`.
3: Unknown or uninitialised column: `germ`.
m1.1 <- ulam(
alist(
day | germ==1 ~ exponential( lambda),
day | germ==0 ~ custom(exponential_lccdf( !Y | lambda)),
lambda <- 1.0 / mu,
log(mu) <- a[temps],
a[temps] ~ normal(0,1)),
data=d,
chains=4,
cores = 4
)
starting worker pid=54244 on localhost:11220 at 11:44:48.270
starting worker pid=54258 on localhost:11220 at 11:44:48.609
starting worker pid=54273 on localhost:11220 at 11:44:49.300
starting worker pid=54288 on localhost:11220 at 11:44:49.607
SAMPLING FOR MODEL '6c7ed290b113c205c3adee6a5a6dddff' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 0.000109 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 1.09 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 1000 [ 0%] (Warmup)
Chain 1: Iteration: 100 / 1000 [ 10%] (Warmup)
Chain 1: Iteration: 200 / 1000 [ 20%] (Warmup)
SAMPLING FOR MODEL '6c7ed290b113c205c3adee6a5a6dddff' NOW (CHAIN 2).
Chain 2:
Chain 2: Gradient evaluation took 7.5e-05 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.75 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2:
Chain 2:
Chain 2: Iteration: 1 / 1000 [ 0%] (Warmup)
Chain 1: Iteration: 300 / 1000 [ 30%] (Warmup)
Chain 2: Iteration: 100 / 1000 [ 10%] (Warmup)
Chain 1: Iteration: 400 / 1000 [ 40%] (Warmup)
Chain 2: Iteration: 200 / 1000 [ 20%] (Warmup)
Chain 1: Iteration: 500 / 1000 [ 50%] (Warmup)
Chain 1: Iteration: 501 / 1000 [ 50%] (Sampling)
SAMPLING FOR MODEL '6c7ed290b113c205c3adee6a5a6dddff' NOW (CHAIN 3).
Chain 3:
Chain 3: Gradient evaluation took 7.4e-05 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.74 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3:
Chain 3:
Chain 3: Iteration: 1 / 1000 [ 0%] (Warmup)
Chain 2: Iteration: 300 / 1000 [ 30%] (Warmup)
Chain 1: Iteration: 600 / 1000 [ 60%] (Sampling)
Chain 2: Iteration: 400 / 1000 [ 40%] (Warmup)
Chain 3: Chain 1: Iteration: 700 / 1000 [ 70%] (Sampling)
Iteration: 100 / 1000 [ 10%] (Warmup)
Chain 2: Iteration: 500 / 1000 [ 50%] (Warmup)
Chain 2: Iteration: 501 / 1000 [ 50%] (Sampling)
Chain 1: Iteration: 800 / 1000 [ 80%] (Sampling)
Chain 3: Iteration: 200 / 1000 [ 20%] (Warmup)
Chain 2: Iteration: 600 / 1000 [ 60%] (Sampling)
Chain 1: Iteration: 900 / 1000 [ 90%] (Sampling)
Chain 2: Iteration: 700 / 1000 [ 70%] (Sampling)
Chain 3: Iteration: 300 / 1000 [ 30%] (Warmup)
Chain 1: Iteration: 1000 / 1000 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 0.294589 seconds (Warm-up)
Chain 1: 0.230801 seconds (Sampling)
Chain 1: 0.52539 seconds (Total)
Chain 1:
Chain 2: Iteration: 800 / 1000 [ 80%] (Sampling)
Chain 3: Iteration: 400 / 1000 [ 40%] (Warmup)
SAMPLING FOR MODEL '6c7ed290b113c205c3adee6a5a6dddff' NOW (CHAIN 4).
Chain 4:
Chain 4: Gradient evaluation took 6.6e-05 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.66 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4:
Chain 4:
Chain 4: Iteration: 1 / 1000 [ 0%] (Warmup)
Chain 2: Iteration: 900 / 1000 [ 90%] (Sampling)
Chain 3: Iteration: 500 / 1000 [ 50%] (Warmup)
Chain 3: Iteration: 501 / 1000 [ 50%] (Sampling)
Chain 2: Iteration: 1000 / 1000 [100%] (Sampling)
Chain 2:
Chain 2: Elapsed Time: 0.280862 seconds (Warm-up)
Chain 2: 0.212282 seconds (Sampling)
Chain 2: 0.493144 seconds (Total)
Chain 2:
Chain 3: Iteration: 600 / 1000 [ 60%] (Sampling)
Chain 4: Iteration: 100 / 1000 [ 10%] (Warmup)
Chain 3: Iteration: 700 / 1000 [ 70%] (Sampling)
Chain 4: Iteration: 200 / 1000 [ 20%] (Warmup)
Chain 3: Iteration: 800 / 1000 [ 80%] (Sampling)
Chain 3: Iteration: 900 / 1000 [ 90%] (Sampling)
Chain 4: Iteration: 300 / 1000 [ 30%] (Warmup)
Chain 3: Iteration: 1000 / 1000 [100%] (Sampling)
Chain 3:
Chain 3: Elapsed Time: 0.310661 seconds (Warm-up)
Chain 3: 0.219301 seconds (Sampling)
Chain 3: 0.529962 seconds (Total)
Chain 3:
Chain 4: Iteration: 400 / 1000 [ 40%] (Warmup)
Chain 4: Iteration: 500 / 1000 [ 50%] (Warmup)
Chain 4: Iteration: 501 / 1000 [ 50%] (Sampling)
Chain 4: Iteration: 600 / 1000 [ 60%] (Sampling)
Chain 4: Iteration: 700 / 1000 [ 70%] (Sampling)
Chain 4: Iteration: 800 / 1000 [ 80%] (Sampling)
Chain 4: Iteration: 900 / 1000 [ 90%] (Sampling)
Chain 4: Iteration: 1000 / 1000 [100%] (Sampling)
Chain 4:
Chain 4: Elapsed Time: 0.323819 seconds (Warm-up)
Chain 4: 0.215416 seconds (Sampling)
Chain 4: 0.539235 seconds (Total)
Chain 4:
precis(m1.1, depth = 2)
Warning messages:
1: Unknown or uninitialised column: `germ`.
2: Unknown or uninitialised column: `germ`.
posterior
preddata <- expand_grid(temps=1:8, realtemps=sort(unique(germ.stdi$temps)), day=1:28)
Warning message:
Unknown or uninitialised column: `germ`.
pred <- link(m1.1, data = preddata)
str(pred)
List of 2
$ lambda: num [1:2000, 1:1792] 0.00439 0.00507 0.00342 0.00459 0.00509 ...
$ mu : num [1:2000, 1:1792] 228 197 293 218 197 ...
mu is average days to germ, lambda is rate parameter. neither change over time, of course, so having days doesn’t really make sense.
preddata$mu <- apply(pred$mu,2,mean)
Warning message:
Unknown or uninitialised column: `germ`.
preddata$low <- apply(pred$mu,2,HPDI)[1,]
preddata$high <- apply(pred$mu, 2, HPDI)[2,]
Single temp. Don’t need day posterior
preddata <- expand_grid(temps=3)
Warning message:
Unknown or uninitialised column: `germ`.
pred <- link(m1.1, data = preddata)
str(pred)
List of 2
$ lambda: num [1:2000, 1] 0.0176 0.0212 0.0178 0.0241 0.0221 ...
$ mu : num [1:2000, 1] 56.9 47.3 56 41.5 45.3 ...
how to convert to probs? use pexp.
predprobs <- pexp(1:28,rate=pred$lambda[1])
Warning messages:
1: Unknown or uninitialised column: `germ`.
2: Unknown or uninitialised column: `germ`.
plot(x=1:28,y=predprobs, type="l")
even thought it isn’t using day, including day in the prediction data frame will help me keep data in the correct format. Maybe.
preddata <- expand_grid(temps=1:8, day=1:28)
Warning message:
Unknown or uninitialised column: `germ`.
pred <- link(m1.1, data = preddata)
str(pred)
List of 2
$ lambda: num [1:2000, 1:224] 0.00439 0.00507 0.00342 0.00459 0.00509 ...
$ mu : num [1:2000, 1:224] 228 197 293 218 197 ...
predresults <- preddata %>%
mutate(lambda=as.list(as.data.frame(pred$lambda)))
predresults
predresults <- predresults %>%
Warning messages:
1: Unknown or uninitialised column: `germ`.
2: Unknown or uninitialised column: `germ`.
mutate(probs=map2(day, lambda, ~ pexp(.x, .y)),
mu=map_dbl(probs, mean),
lower=map_dbl(probs, ~ HPDI(.)[1] %>% unlist()),
upper=map_dbl(probs, ~ HPDI(.)[2]) %>% unlist())
predresults
predresults %>% select(-lambda, -probs) %>%
Warning messages:
1: Unknown or uninitialised column: `germ`.
2: Unknown or uninitialised column: `germ`.
3: Unknown or uninitialised column: `germ`.
4: Unknown or uninitialised column: `germ`.
5: Unknown or uninitialised column: `germ`.
6: Unknown or uninitialised column: `germ`.
7: Unknown or uninitialised column: `germ`.
8: Unknown or uninitialised column: `germ`.
9: Unknown or uninitialised column: `germ`.
mutate(temps=factor(temps, labels=as.character(sort(unique(germ.stdi$temps))))) %>%
ggplot(aes(x=day,y=mu,color=temps,group=temps)) +
geom_line()
Add realdata:
stdi.plot <- germ %>% filter(pops=="STDI") %>%
Warning messages:
1: Unknown or uninitialised column: `germ`.
2: Unknown or uninitialised column: `germ`.
3: Unknown or uninitialised column: `germ`.
select(day, temps, cumulative_germ, total_seeds) %>%
mutate(temps=as.factor(temps),
prop.germ=cumulative_germ/total_seeds)
predresults %>% select(-lambda, -probs) %>%
mutate(temps=factor(temps, labels=as.character(sort(unique(germ.stdi$temps))))) %>%
ggplot(aes(x=day,y=mu,color=temps,group=temps)) +
geom_line() +
geom_point(aes(y=prop.germ), data=stdi.plot)
Poor!